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1.  INTRODUCTION 


The  work  presented  considers  the  adaptation  of  potential  theory  and  the  boundary  integral 
method  for  the  development  of  computer  codes  to  simulate  underwater  explosion  bubble 
dynamics.  In  order  to  understand  the  development  of  numerical  methods  like  this,  it  is 
important  to  review  some  of  the  pioneering  efforts  which  have  lead  to  its  use.  One  of  the 
early  studies  in  the  area  of  bubble  dynamics  was  published  by  Benjamin  and  Ellis  (1966). 
Much  of  the  work  that  followed  builds  from  this  fundamental  study.  Benjamin  and  Ellis’  work 
provided  a  broad  examination  of  the  subject  matter  referring  to  the  well-known  implosion 
mechanism  first  recognized  by  Lord  Rayleigh  (1917).  Their  study  included  a  review  of  the 
contributions  of  such  notable  authors  as  Taylor  (1950),  Taylor  and  Davies  (1950),  Cole  (1948), 
Lamb  (1945),  Knapp  and  Hollander  (1948),  Knapp  (1958a,  1958b),  and  Hammitt  (1963).  The 
effort  also  provided  insight  into  the  physics  of  bubble  jetting.  In  doing  so,  the  mechanisms 
leading  to  bubble  collapse  and  jetting  were  examined  from  both  a  physical  and  mathematical 
standpoint.  In  their  presentation,  the  authors  considered  the  importance  of  the  Kelvin  impulse. 
The  bubble's  impulse  was  used  to  explain  the  nature  of  bubble  collapse,  jetting,  and  the 
formation  of  toroidal  bubbles.  They  explained,  "As  a  general  interpretation  of  the  phenomenon 
of  jet  penetration  through  collapsing  cavities,  one  can  simply  say  that  the  liquid  is  finding  the 
only  possible  way  to  preserve  its  impulse  as  the  cavity  size  decreases."  After  a  thorough 
overview  of  the  bubble  pulsation  and  jetting  phenomenon,  the  authors  proceeded  to 
investigate  the  phenomenon  through  a  series  of  experiments.  The  primary  concern  of  their 
study  was  with  the  small-scale  hydrodynamic  processes  that  lead  to  cavitation  damage. 
Therefore,  the  experimental  investigation  was  initiated  with  the  aim  of  examining  the 
phenomena  of  bubble  collapse  and  jetting  near  rigid  bodies.  In  the  experiments,  a  free-fall 
apparatus  was  used  which  created  a  gravity-free  environment  to  observe  bubble  collapse. 
Additionally,  by  holding  the  experimental  apparatus  fixed,  the  effects  of  gravity  on  bubble 
collapse  were  observed.  Using  this  device,  the  effects  of  a  nearby  boundary  and  gravity  on 
bubble  jetting  were  observed  separately.  The  results  provided  conclusive  proof  of  bubble 
jetting  in  the  free  field,  due  to  gravity,  as  well  as  bubble  jetting  influenced  by  a  nearby  rigid 
wall. 

Studying  the  same  phenomena,  Piesset  and  Chapman  (1971)  developed  numerical 
techniques  designed  especially  for  the  calculation  of  bubble  collapse  and  jetting  near  rigid 
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toundaries.  In  this  model,  the  effects  of  viscosity  and  compressibility  of  the  fluid  were 
legiected  in  the  basic  conservation  laws  for  mass  and  momentum.  A  number  of  other 
amplifying  assumptions  were  made  in  their  model.  In  summary,  they  are  as  follows: 

(1)  The  liquid  is  incompressible. 

(2)  The  flow  is  inviscid. 

(3)  The  vapor  pressure  inside  the  bubble  is  uniform  and  constant. 

(4)  The  ambient  fluid  pressure  remains  constant  with  time. 

(5)  The  bubble  contains  no  permanent  gas. 

(6)  Surface  tension  effects  are  negligible. 

(7)  The  effects  of  gravity  are  not  included. 

(8)  The  bubble  is  initially  spherical  and  remains  axisymmetric  thereafter. 

As  pointed  out  in  their  paper,  only  the  first  three  assumptions  are  essential  to  the  method 
developed.  Correspondingly,  the  same  three  initial  assumptions  are  adopted  in  the  current 
study.  Under  the  assumptions  used  in  their  model,  the  velocity  is  equal  to  the  gradient  of  a 
potential  which  satisfies  Laplace's  equation,  and  the  pressure  is  given  by  Bernoulli’s  equation. 
The  analysis  was  carried  out  based  on  these  equations  and  a  finite  difference  integration 
scheme.  In  their  investigation,  Plesset  and  Chapman  simulated  two  specific  cases  of  an 
initially  spherical  bubble  collapsing  near  a  plane  solid  wall.  Their  calculation  predicted  bubble 
jet  speeds  sufficient  to  explain  cavitation  damage. 

In  a  related  numerical  study,  Chapman  and  Plesset  (1972)  extended  their  original 
calculations  to  include  vapor  bubbles  in  an  infinite  liquid  initially  perturbed  from  the  spherical 
shape.  In  this  study,  two  cases  of  bubble  collapse  were  simulated,  namely,  an  oblate  ellipsoid 
and  a  prolate  ellipsoid.  The  calculations  showed  the  formation  of  a  high-speed  bubble  jet 
during  the  late  stages  of  bubble  collapse.  In  both  calculations,  the  bubble’s  initially  elliptical 
shape  was  transformed  into  a  dumbbell  shape  in  the  opposite  direction  of  the  original 
elongation  during  the  bubble’s  collapse.  The  results  also  showed  the  important  point— that  a 
very  small,  initial  deformation  from  a  spherical  shrpe  limits  the  applicability  of  the  assumption 
of  spherical  symmetry.  One  final  important  point  was  that  the  analysis  indicated  that  the  jet 
speed  at  the  time  of  impact  is  smaller  than  had  been  found  using  linearized  theories. 
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Another  theoretical  study  important  to  bubble  dynamics  was  presented  by  Blake  and 
Gibson  (1981).  In  their  numerical  approach,  they  proceeded  under  the  same  fundamental 
assumptions  of  Plesset  and  Chapman  (1971)  using  incompressible,  irrotational  flow  and  a 
uniform,  internal  bubble  pressure.  The  solution  centered  around  the  use  of  a  number  of  ring 
distributions  of  singularities  to  calculate  the  potential.  The  new  method  was  physically  more 
appropriate  than  previous  methods  based  on  line  distributions  (Bevir  and  Fielding  1974).  In 
their  paper,  Blake  and  Gibson  investigated  the  damaging  effects  of  cavitation  bubbles  and 
methods  to  reduce  the  possibility  of  cavitation  damage.  As  pointed  out  in  their  work  and  from 
the  law  of  Bjerknes,  oscillating  bubbles  tend  to  migrate  toward  rigid  boundaries  and  away  from 
a  free  surface.  Furthermore,  the  migration  and  jetting  of  bubbles  near  flexible  boundaries  are 
not  fully  understood.  Their  investigation  into  flexible  surfaces  concentrated  around  both 
numerical  and  experimental  studies  using  spark-generated  bubbles  near  a  free  surface — a 
free  surface  being  the  limit  of  flexible  surfaces.  The  study  revealed  that  the  approximate 
integral-equation  approach  presented  was  adequate  to  simulate  the  growth  and  collapse  of 
vapor  bubbles  at  least  one  maximum  bubble  radius  from  the  free  surface.  They  additionally 
observed  that  vapor  bubbles  in  the  absence  of  gravity  can  be  generated  within  half  a 
maximum  radius  of  the  free  surface  without  venting.  Their  work  provided  some  excellent 
photography  and  insight  into  the  behavior  of  bubbles  near  a  free  surface. 

An  alternate  numerical  approach  which  was  published  around  the  same  time  as  the  Blake 
and  Gibson  (1981)  contribution  and  was  written  by  Guerri,  Lucca,  and  Prosperetti  (1980). 

Their  method  entailed  the  use  of  the  now-familiar  boundary  integral  numerical  method  for  the 
simulation  of  nonspherical  cavitation  bubbles  in  invisdd,  incompressible  liquids.  Their  work 
offered  an  attractive  alternative  to  the  finite  difference  (Plesset  and  Chapman  1971), 
perturbation  (Chapman  and  Plesset  1972;  Chahine  1977),  and  singularity  (Blake  and  Gibson 
1981)  methods  that  had  been  previously  presented.  The  crucial  step  in  the  procedure 
involves  the  solution  of  the  potential  problem.  The  potential  solution  is  found  using  the 
boundary  integral  or  boundary  element  method  which  is  based  on  an  exact  relation,  namely, 
Green's  identity.  One  of  the  more  attractive  features  of  their  approach  is  the  method’s 
efficiency.  The  efficiency  comes  from  the  need  of  only  surface  values  for  the  velocity  potential 
and  its  first  derivative.  Therefore,  the  method  avoids  the  problem  of  solving  Laplace’s 
equation  over  the  entire  domain  occupied  by  the  liquid.  Another  appealing  feature  of  the 
boundary  integral  method  is  its  excellent  performance  in  calculating  the  unstable  collapse 


3 


phase  of  the  bubble  without  smoothing  techniques.  In  order  to  show  the  reliability  of  the 
computational  model,  one  of  the  cases  originally  studied  by  Plesset  and  Chapman  was 
considered  for  comparison.  In  this  comparison,  the  velocities  of  the  bubble  jet  tip  agreed 
favorably  with  the  previous  calculations  of  Plesset  and  Chapman.  The  results  presented 
several  examples  of  the  method’s  capability  to  predict  bubble  collapse  and  jetting  in 
axisymmetric  configurations.  They  included  the  collapse  of  a  bubble  near  a  rigid  wall,  as 
previously  studied  by  Plesset  and  Chapman,  and  the  calculation  of  bubble  collapse  for  three 
bubbles  with  colinear  centers.  Their  method  was  crucial  for  future  multidimensional  bubble 
codes. 

Another  major  contribution  in  the  field  of  bubble  dynamics  was  presented  by  Chahine, 
Perdue,  and  Tucker  (1988).  Their  work  introduced,  for  the  first  time,  an  economical  approach 
to  the  numerical  study  of  three-dimensional  (3-D)  bubble  dynamics.  The  essence  of  the 
method  was  based  on  the  boundary  integral  approach.  Nonetheless,  their  study  revealed  a 
number  of  limiting  factors  in  the  boundary  integral  approach’s  applicability  for  3-D  bubble 
dynamics,  as  well  as  a  number  of  solutions  to  unforeseen  problems  in  the  method.  Their 
work  required  a  number  of  contributions  in  the  process  of  developing  numerical  methods  to 
simulate  underwater  explosion  bubble  dynamics.  One  of  the  hurdles  that  Chahine  and  Perdue 
needed  to  overcome  was  the  calculation  of  the  boundary  integral  method’s  surface  integrals  in 
3-D  configurations.  Originally,  the  numerical  calculation  of  these  integrals  was  thought  too 
formidable  for  simple  numerical  approximation.  Therefore,  Perdue  (1988)  underwent  the 
laborious  task  of  calculating  the  integrals  out  exactly.  Chahine  and  Perdue  also  devised  some 
innovative  methods  for  the  calculation  of  the  tangential  components  of  velocity  on  the  surface 
of  the  bubble.  In  all  these  approaches,  a  number  of  unforeseen  problems  with  3-D  boundary 
integral  techniques  arose.  They  included  geometric  perturbations  caused  by  the  choice  of  the 
initial  grid,  solid  angle  calculations,  initial  geometries,  and  numerical  instabilities,  to  mention 
only  a  few.  In  Chahine  and  Perdue’s  report  on  the  approach,  they  focused  on  explosion 
bubbles  as  their  topic  of  interest  and  discussed  the  application  of  their  work  in  calculating 
bubble  collapse  near  submerged  bodies.  They  also  presented  a  number  of  calculations  with 
rigid  plates  at  various  angles  of  inclination,  a  calculation  near  a  submerged  body,  and  some 
axisymmetric  calculations.  Their  work  represented  a  significant  contribution  in  the  field  of 
computational  bubble  dynamics.  Their  work  additionally  offered  a  springboard  into  numerous 
unanswered  questions  regarding  the  boundary  integral  method  and  3-D  computational  models. 
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2.  MATHEMATICAL  DEVELOPMENT 


The  first  step  in  the  formulation  of  this  bubble  dynamics  model  was  to  put  the  basic  laws 
of  nature  governing  the  phenomenon  into  a  suitable  mathematical  form.  The  physical 
conditions  imposed  on  the  model  are  depicted  graphically  in  Figure  1 .  The  fluid  occupies  a 
domain  ft,  bounded  by  the  bubble  surface,  Sb,  solid  boundaries,  Sr,  and  a  surface  at  infinity, 
S_.  The  fluid  contained  in  the  domain  ft  is  considered  incompressible,  and  inviscid  flow  is 
also  assumed.  Furthermore,  the  gas  in  the  bubble  is  assumed  adiabatic.  In  order  to  justify 
these  assumptions,  the  physical  nature  of  an  underwater  explosion  bubble  must  be 
considered.  In  the  case  of  underwater  explosive  detonations,  a  chemical  reaction  converts 
the  original  material  very  quickly  into  a  gas  at  a  very  high  temperature  and  pressure.  The 
temperature  of  the  explosive  products  directly  after  detonation  can  be  of  the  order  of  3,000°  0. 
with  pressures  on  the  order  of  50,000  mpa  (Cole  1946).  The  result  of  the  initial  detonation  is 
a  compression  or  shock  wave  being  emitted  into  the  fluid  field,  followed  by  the  dynamic 
expansion  and  contraction  of  the  gaseous  products.  After  the  release  of  the  shock  wave, 
which  is  an  early  time  phenomenon,  the  speed  of  the  bubble's  surface  remains  at  least  an 
order  of  magnitude  smaller  than  the  speed  of  sound  in  water.  Therefore,  the  imposition  of  an 
incompressible  condition  on  the  fluid  is  deemed  valid  for  this  flow  after  the  first  few 
microseconds. 

Using  the  bubble  radius  in  the  expansion  or  contraction  phase  of  the  bubble  pulse  as  a 
characteristic  length,  a  Reynolds  number  can  be  calculated.  This  Reynolds  number  is  found 
to  be  high  through  most  of  the  bubble  growth  and  collapse.  Recalling  that  the  Reynolds 
number  is  the  ratio  of  inertial  to  viscous  forces,  it  is  clear  that  neglecting  viscous  terms  in  the 
conservation  of  momentum  equation  detracts  very  little  from  the  accuracy  of  the  solution. 
Additionally,  under  the  assumptions  that  the  fluid  flow  is  inviscid  and  irrotational,  the  velocity 
field  can  be  found  from  Laplace’s  equation.  In  the  case  of  a  cavitation  bubble,  similar 
assumptions  were  shown  to  be  valid  by  Guerri,  Lucca,  and  Prosperetti  (1980)  and  others 
(Blake,  Taib,  and  Doherty  1986;  Chahine  1982;  Oguz  and  Prosperetti  1989).  In  summary,  the 
gradient  of  the  potential  normal  to  a  rigid  body  is  zero,  the  potential  can  be  determined  on  the 
bubble  surface  initially,  and  the  velocity  and  potential  at  infinity  vanish. 
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Figure  1 .  Definition  of  the  Physical  Problem. 


The  physical  problem  and  the  assumptions  described  previously  can  be  stated  in 
mathematical  terms.  The  resulting  mathematical  formulation  consists  of  two  coupled 
equations,  one  derived  from  conservation  of  momentum  and  the  other  from  the  assumptions 
of  irrotational  flow  and  conservation  of  mass.  Under  the  assumptions  stated,  these 
fundamental  equations  take  the  form  of  Bernoulli’s  equation  and  Laplace’s  equations.  The 
solution  of  these  two  coupled  equations  is  the  basis  of  the  computational  model. 


The  first  assumption  is  that  the  liquid  is  incompressible  so  that  the  divergence  of  the 
velocity  vector  is  zero.  This  result  derives  from  the  conservation  of  mass  equation  and  is  as 
follows: 

V-u-0.  (!) 

For  irrotational  flow,  the  curl  of  the  velocity  field  is  zero, 

V  X  u  -  0.  (2) 

The  previous  expression  is  only  satisfied  when  the  velocity  u  is  a  function  of  the  potential 
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u  *  . 


(3) 


Combining  Equations  1  and  3  yields  Laplace’s  equation  to  be  satisfied  in  the  domain  a 
occupied  by  the  fluid, 

V24>  -  0.  W 

The  boundary  conditions  under  which  the  equation  must  be  solved  are  as  follows: 

on  Sb ,  (5) 

dt 


4»  ->  0  |  r  |  , 


(6) 


i!  -  V$  -n  *  0  on  Sf  ,  (7) 

dn 


where  r  is  the  position  vector  of  a  generic  field  point  The  pressure  in  the  liquid  is  given  by 
Bernoulli’s  equation  in  the  following: 


3* 

at 


(8) 


where  p  is  the  liquid  density,  g  is  acceleration  of  gravity,  and  z  is  the  vertical  component  of 
the  position  vector  r .  In  Equation  8,  P.  denotes  the  pressure  at  a  large  distance  from  the 
bubble  on  the  plane  z=0.  Insertion  of  the  total  derivative  cfy/dt  =  d$/9t  +  u  •  V$  =  a$/at  +  u  2 
into  Equation  8  gives  the  following: 


d$  u*  P 
—I.  *  — *  ♦  _ 

dt  2  p 


(9) 


The  dynamic  boundary  condition  on  the  surface  of  the  bubble  requires  the  liquid  pressure  to 
equal  the  internal  pressure  of  the  gas  Pg.  With  this,  integration  in  time  of  the  expression  in 
Equation  9  gives  a  boundary  condition  on  the  bubble's  surface  for  the  potential.  Furthermore, 
the  tangential  derivatives  of  the  potential  on  the  bubble's  surface  can  be  determined  from 
Equation  9.  The  remaining  quantities  to  be  determined  are  the  gradient  of  the  potential 
normal  to  the  surface  on  the  bubble  and  the  internal  gas  pressure.  Thus,  Equations  4,  9,  and 
the  boundary  conditions  stated  in  Equations  5-7  form  the  basis  for  the  current  computational 
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model.  The  sections  to  follow  will  describe  the  mathematical  model  required  for  the  solution 
of  these  two  coupled  equations  under  the  boundary  conditions  prescribed. 


One  of  the  terms  requiring  further  consideration  in  the  conservation  of  momentum  equation 
is  the  pressure  at  points  on  the  surface  of  the  bubble.  Here,  the  pressure  must  equal  the  sum 
of  the  internal  gas  pressure  and  a  surface  tension  component.  For  explosion  bubbles,  the 
surface  tension  component  has  been  shown  to  be  of  little  importance  for  most  cases  of 
explosion  bubble  dynamics  (Wilkerson  1989, 1990).  However,  for  cavitation  and  gas  bubbles, 
this  component  can  be  of  importance.  The  relation  for  the  pressure  on  the  surface  of  an 
explosion  bubble  can  be  written  as  follows: 

P  -  P9  .  <10) 


where  PB  is  the  internal  gas  pressure  of  the  bubble.  This  pressure,  PB,  must  be  determined 
from  an  equation  of  state.  After  the  initial  gas  pressure  is  found,  the  internal  gas  pressure  can 
be  updated  from  the  ratio  of  volumes  as  follows: 


P 


8 


go 


\7 

> 


(11) 


where  Vc  is  the  current  volume  V,  is  the  initial  volume,  PB0  is  the  initial  internal  gas  pressure, 
and  y  is  the  adiabatic  constant  which  depends  on  the  explosive  type.  The  adiabatic  law  is 
approximate  because  the  thermal  penetration  length  in  the  bubble  is  much  smaller  than  its 
radius.  The  adiabatic  constant  y  has  been  empirically  extracted  from  underwater  explosive 

testing. 


Effective  initial  conditions  are  needed  to  determine  P^  and  the  initial  radius.  As  a  first 
guess  for  the  initial  internal  bubble  gas  pressure  Pgo,  the  mass  of  the  explosive  charge  can  be 
converted  directly  into  a  high  pressure  gas,  and  an  adiabatic  pressure  volume  relation  can  be 
used  thereafter.  However,  due  to  the  finite  time  required  for  the  explosive  to  convert  into  a 
gas,  the  effects  of  compressibility  in  the  fluid,  and  the  small  expansion  of  the  gas  products 
during  the  chemical  reaction,  this  assumption  results  in  a  poor  initial  condition  for  the  pressure 
and  the  initial  bubble  radius.  The  resulting  calculation  for  predicting  bubble  radius  and  bubble 
periods  has  unacceptably  high  errors  when  compared  to  experimental  results. 
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Recently,  Chahine  (1989)  has  shown  that  with  experimental  knowledge  of  minimum  and 
maximum  bubble  radii  for  an  explosive  type,  the  Rayleigh-Plesset  equation  can  be  integrated, 
yielding  more  consistent  initial  conditions  for  the  internal  pressure  at  a  given  radius.  The 
relation  discussed  in  Chahine  (1989)  is  based  on  the  Rayleigh-Plesset  equation  or  spherical 
bubbles  and  some  empirical  data  from  explosion  bubble  experimentation.  The  principles  of 
the  expression  are  given  here.  The  Rayleigh-Plesset  equation  is  as  follows: 


p[Rft 


(12) 


where  R  is  the  radius  and  P,  is  the  ambient  pressure  far  away  from  the  bubble  at  the  initial 
depth  of  its  center.  Here,  the  dots  denote  derivatives  with  respect  to  time.  When  the  ambient 
pressure  P.  is  constant.  Equation  12  is  integrate  when  written  in  the  following  form  (Plesset 
1949): 

-  P,(Ro9r)R*‘*r  -  P.R2.  (13) 

Integrating  Equation  13  and  allowing  the  current  radius  R  to  equal  the  maximum  bubble 
radius,  R0  to  equal  the  minimum  or  initial  bubble  radius,  and  the  corresponding  velocities  to 
equal  zero,  yields  the  following  relation  for  the  initial  bubble  pressure: 


£  f_5L(R3P 

2  1  dR ' 


R2) 


P 


go 


P.(Y-1)(1-cr») 
1  » 


(14) 


where  a  =  R^R^  is  the  ratio  of  minimum  to  maximum  bubble  radii.  For  a  more  detailed 
derivation  of  Equation  14,  see  Chahine  (1 989).  Data  for  the  ratio  a  can  be  found  for  a 
number  of  explosive  types  in  experimental  studies  (Snay,  Goertner,  and  Price  1952).  This 
relation  applies  at  large  depths  where  an  explosion  bubble  remains  spherical.  Since  bubbles 
are  approximately  spherical  during  their  growth  in  shallow  water  as  well,  this  relation  applies 
approximately  by  adjusting  P,.  Therefore,  the  remaining  problem  is  how  to  obtain  the  initial 
radius.  The  initial  radius  can  be  determined  from  one  of  two  methods.  One  method  involves 
an  empirical  relation  for  the  pressure  P^  based  on  the  weight,  explosive  type,  and  volume  of 
the  charge.  Once  such  accepted  pressure  relation  for  TNT  is  (Cole  1948)  as  follows: 
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(15) 


where  W  is  the  mass  of  TNT  in  grams,  V  is  the  volume  in  cm3,  P#0  is  the  pressure  in  kilobars, 
and  y  is  the  adiabatic  constant  for  TNT,  which  is  approximately  1 .25.  This  relation  is  valid  for 
P  <  4,500  lb/in2  pressure.  The  previous  relation  has  been  confirmed  for  TNT  underwater 
explosive  by  Wilkerson  (1989, 1990). 


Another  method  of  determining  the  initial  pressure  is  based  on  an  empirical  relation  for  the 
maximum  radius  and  the  previous  expression  relating  the  minimum  to  maximum  bubble 
radii  a.  This  approach  requires  an  estimation  of  Rt71„.  This  can  be  determined  on  the  basis  of 
empirical  rules  for  explosion  bubbles.  This  rule  came  from  a  dimensional  analysis  argument 
which  related  explosive  energy  to  bubble  energy.  While  the  exact  roots  of  the  expressions 
given  here  are  difficult  to  trace,  the  formula  for  maximum  bubble  radius  of  an  explosive  charge 
can  be  deduced  from  the  following  argument.  The  energy  of  the  explosive  charge  E  is 
proportional  to  the  product  of  the  maximum  volume  of  water  displaced  at  a  given  hydrostatic 
pressure  times  a  constant  K  as  follows: 

E  =  KPV.  (16) 


The  energy  of  the  charge  is  proportional  to  the  weight  of  the  charge  W,  the  hydrostatic 
pressure  is  proportional  to  the  depth  Z  plus  the  ambient  pressure  at  the  fluid’s  surface  (e.g., 
D  =  Z  +  33  where  all  depths  are  in  feet),  and  the  maximum  volume  is  proportional  to  the 
maximum  radius  of  the  bubble  cubed  R3^.  (At  a  depth  of  33  ft,  the  hydrostatic  head  equals 
approximately  1  atmosphere  of  pressure.)  From  these  relations,  an  expression  for  the 
maximum  bubble  radius  can  be  written  as  follows: 


(17) 


where  J  is  a  constant  to  be  determined  from  experimental  results.  The  constant  J  has  been 
determined  for  a  number  of  explosive  types  at  deep  depths  from  underwater  explosion 
experimentation  (Swift  and  Decius  1950). 
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The  remaining  quantities  to  be  prescribed  for  the  boundary  integral  method  are  the  initial 
values  of  the  velocity  and  potential.  The  velocity  is  assumed  to  vanish,  u  =  0  at  t  =  0.  Since 
the  physically  significant  quantity  is  u  =  V$,  4>  can  be  taken  to  vanish  initially. 


To  simplify  the  validation  of  the  codes,  it  is  convenient  to  use  scaling  factors,  making  the 
results  nondimensional.  The  parameters  chosen  are  consistent  with  the  thinking  used  to 
determine  the  initial  conditions  for  the  model.  P,  =  hydrostatic  plus  atmospheric  pressure  at 
the  assumed  depth  of  the  initial  bubble  (e.g.,  an  explosive  at  a  33-ft  depth  would  have  a  P,  of 
approximately  two  atmospheres);  p  is  the  density  of  the  fluid;  and  is  given  by  Equation  17 
and  corresponds  to  the  observed  maximum  radius  of  an  explosion  bubble  in  the  free  field  (the 
empirical  rule  just  derived  for  maximum  bubble  radius  in  the  absence  of  gravity).  Using  these 
three  parameters,  scale  factors  for  time  T,,  length  L,,  and  mass  M,  were  defined.  They  are  as 
follows: 

TS  *  Rm«  l/p/P.  • 

L.  -  Rm«  .  <18> 

Mi  *  • 


Dimensionless  pressures  are  defined  by  P/P,.  In  terms  of  these  fundamental  quantities,  the 
mathematical  model  may  be  written  as  follows: 


V2  *  0  . 


d4> 

"dt 


(19) 


where  all  quantities  are  dimensionless,  and  the  Froude  number  is  defined  by  Fr  =  P^pgL,  = 
d/Rm...  where  d  is  the  depth  of  the  charge.  It  is  seen  from  Equation  19  that  for  small  values  of 
F„  the  effect  of  gravity  is  important,  while  it  is  negligible  for  Fr »  1 .  in  particular,  the  limit  in 
which  the  Rayleigh-Plesset  equation  is  appropriate  is  obtained  for  Fr  «>. 
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The  boundary  integral  formulation  for  the  solution  of  Laplace’s  equation  is  based  on  the 
use  of  Green’s  second  identity  which  reduces  a  volume  integral  to  a  surface  integral.  The 
remaining  surface  integral  is  discretized,  yielding  a  linear  algebraic  system  for  the  gradient  of 
the  potential  normal  to  the  boundaries.  The  starting  point  of  the  boundary  integral  method  is 
Green's  second  identity  as  follows: 

XfoV’y-yV2*)^- -y|£)dS.  (20) 


which  is  valid  for  any  pair  of  sufficiently  smooth  functions  4>  and  y.  In  the  present  approach,  $ 
is  taken  to  be  the  velocity  potential  discussed  in  Section  2.2,  and  y  is  taken  to  be  the  free 
space  Green’s  function  for  Laplace’s  equation  as  follows: 


V 


1 

l£  -  R  |* 


(21) 


Physically,  this  represents  the  field  at  position  JR,  generated  by  a  source  of  intensity  4n  placed 
at  r .  This  function  satisfies  an  equation  of  the  following  form: 


V2y  «  -  X8(r  -  R)  , 


(22) 


where  the  Laplacian  is  with  respect  to  the  variable  R.  Here,  X  =  4tt  if  R  is  internal  to  the 
domain  ft.  When  R  is  located  on  the  boundary  (in  particular,  on  the  bubble’s  surface  S),  X 
equals  the  solid  angle  under  which  the  domain  ft  is  viewed  from  R.  It  is  important  to  note 
that  X  must  be  evaluated  for  the  actual  discretized  shape  of  the  boundary  rather  than  for  the 
exact  one;  otherwise,  incorrect  numerical  results  would  be  found.  Substitution  of  Equations  21 
and  22  into  Equation  20  yields  the  following: 


X*(R) 


1 

iL-m 


d<t>(r) 

dnr 


1 

il-ei 


dS  • 


(23) 


The  model  updating  in  this  formulation  involves  the  time  integration  of  Equation  9  and  the 
tracking  of  the  bubble's  surface  with  time.  In  order  to  follow  the  bubble’s  shape,  an  accurate 
estimation  of  the  bubble’s  surface  velocity,  based  on  the  tangential  and  normal  velocity,  must 
be  determined.  Inaccuracies  in  the  estimation  of  the  bubble's  velocity  must  be  determined. 
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Inaccuracies  in  the  estimation  of  the  bubble's  velocity  will  add  to  instabilities  in  the  numerical 
solution.  As  a  first  attempt,  a  first-order  finite  difference  scheme  was  used.  However,  the 
resulting  integration  required  a  small  time  step  and  caused  the  overall  computation  to  be  too 
long.  Therefore,  a  combination  of  an  explicit  and  implicit  finite  difference  scheme  was  used  to 
help  stabilize  the  bubble  code’s  time  integration  and  model  updating.  The  combination  of  the 
two  expressions  resulted  in  a  method  accurate  to  the  second  order  and  far  more  stable 
numerically.  The  method  used  has  the  following  form: 

xM  =x,  +  (x,+xw)At/2  ,  (24) 

for  each  of  the  Cartesian  coordinates,  and 


♦u  =  +  (0.  +  (25) 

for  updating  the  potential.  Other  possible  methods  to  help  stabilize  the  time  integration  also 
exist.  One  such  method  is  given  in  Chahine,  Perdue,  and  Tucker  (1988).  In  this  approach,  a 
term  for  artificial  viscosity  is  included  in  Equation  9.  This  term  accomplishes  the  same  result 
as  the  combination  of  the  finite  difference  schemes.  Another  possibility  would  be  to  use  a 
higher-order  predictor  corrector  method.  This  method  would  significantly  increase  computation 
time,  thereby  defeating  some  of  the  appeal  of  the  current  method. 

For  a  single  degree  of  freedom  system,  namely,  a  spherically  symmetric  system, 

Equation  12  can  be  integrated  directly  for  a  solution.  Such  solutions  have  been  shown  in 
good  agreement  with  experimental  and  other  numerical  studies  (Wilkerson  1989).  For  an 
axisymmetric  configuration,  it  is  necessary  to  break  up  the  bubble's  boundary  into  discrete 
segments.  For  the  current  development,  the  initial  axisymmetric  discretization  consists  of 
equally  spaced,  straight  line  segments  around  the  bubble’s  shape.  This  initial  discretization, 
shown  in  Figure  2,  shows  how  the  bubble's  shape  can  be  approximated  by  n  straight  line 
segments  with  n  +  1  nodes.  The  geometry  is  based  on  a  cylindrical  coordinate  system  (r,  z,0) 
where  8  is  the  counterclockwise  rotation  about  the  z  axis. 

For  the  axisymmetric  formulation,  the  vectors  r  and  R  locating  field  points  on  boundaries 
are  given  by  the  following  relations: 
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Figure  2.  Axisymmetric  Discretization  (Cylindrical  Coordinates). 


r  «  [rcos(0),rsin(0),z]  , 
R  «  [R,0,Z]  . 


(26) 


The  3-D  formulation  offered  a  number  of  alternatives  with  regards  to  the  initial  gridding.  Some 
of  the  possibilities  included  three-noded  triangles,  rectangular  patches,  or  other  shaped 
elements.  For  the  current  numerical  model,  three-noded  triangular  patches  were  used.  The 
first  level  of  discretization  for  the  3-D  model  is  based  on  the  initial  shape  of  an  icosahedron. 
This  shape  consists  of  20  equilateral  triangles  and  12  nodes,  all  of  which  lie  on  the  surface  of 
a  sphere.  This  initial  level  of  discretization  can  be  improved  by  breaking  up  each  of  the 
original  triangles  into  smaller  triangles.  With  each  additional  level  of  refinement,  the 
coarseness  of  the  original  icosahedron  shape  can  be  gradually  improved  to  more  closely 
approximate  a  near-spherical  shape.  The  process  of  level  refinement  is  shown  graphically  in 
Figure  3. 
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Level  1 
20  Elements 
12  Nodes 


Level  2 
80  Elements 
42  Nodes 


Level  3 
180  Elements 
92  Nodes 


Figure  3.  First  Three  Levels  of  Discretization  for  Triangles. 

A  field  point  anywhere  on  the  surface  of  an  element  can  be  defined  by  the  position  vector 
r  as  follows: 


9  3  3 

I  -  £  N,x,i  +  £  N,y,j  +  £  N,z,k  . 

WWW 


(27) 


where  N,  is  the  shape  function  for  a  standard  triangular  element,  x„  y„  z,  are  the  nodal 
coordinates,  and  i,  j,  k  the  Cartesian  coordinate  system’s  unit  vectors.  Since  these 
relationships  apply  to  any  linear  quantity  on  the  surface  of  the  element,  a  similar  method  can 
be  applied  for  interpolating  the  potential  and  gradient  of  the  potential  on  the  surface  in  terms 
of  the  shape  functions  and  the  nodal  values  of  these  quantities.  These  relations  are  as 
follows: 


3 


4  *  £  ni$i  • 

to 


(28) 


15 


(29) 


9<t> 

3n 


\ 

> 


3 


-  £  N, 
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where  N,  is  the  same  shape  function  given  in  Equation  27,  $  and  3<t>/3n  are  the  interpolated 
potential  and  gradient  of  the  potential  anywhere  on  the  surface  of  the  element,  and  <$>,  and 
d$/dn  are  the  nodal  values  of  the  potential  and  gradient  of  the  potential.  The  global 
coordinates  of  the  R  vector  are  as  follows: 

R,  =  Rx,i  ♦  Ryli  +  Rzlk  .  (30) 

where  i  refers  to  the  global  node  number. 

3.  RESULTS 

Comparisons  were  made  with  the  axisymmetric  formulation  for  free-field  bubble  jetting. 
The  present  axisymmetric  code  has  been  shown  to  be  in  good  agreement  with  physical 
phenomena  and  other  numerical  approaches  for  the  calculation  of  bubble  jet  velocities  by 
Wilkerson  (1989).  Further,  direct  comparisons  with  the  axisymmetric  code  of  Oguz  and 
Prosperetti  (1989)  were  also  made. 

The  first  comparison  is  between  the  axisymmetric  code  and  the  first-order  3-D  code.  This 
comparison  uses  the  results  from  a  calculation  for  a  1  -lb  charge  of  TNT  with  a  total 
hydrostatic  head  of  33  ft.  Under  these  conditions,  there  is  a  strong  influence  of  gravity,  and 
the  bubble  should  form  a  jet  firing  the  collapse  phase  of  its  first  pulsation.  Plots  comparing 
velocity  and  displacement  of  the  north  and  south  poles  of  the  bubble  are  provided.  Figure  4 
compares  the  absolute  velocity  of  the  bubble’s  north  pole  from  the  two  calculations.  Similarly, 
Figure  5  compares  the  absolute  velocity  of  *he  bubble’s  south  poles.  It  is  important  to  note 
that  these  calculations  were  made  using  dimensionless  quantities.  Note  that  for  the 
dimensionless  units  used,  the  fundamental  velocity  is  y/p,/p  =  33.0  ft/s.  Therefore,  the 
velocities  from  the  south  pole,  where  the  bubble  jet  is  formed  during  the  late  stages  of 
collapse,  are  quite  large. 
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Figure  4.  Comparison  of  Bubble  North  Pole  Surface  Velocity  as  Computed  by  the 
Axisvmmetric  and  3-D  Bubble  Codes. 


Figure  5.  Comparisons  of  Bubble  South  Pole  Surface  Velocity  as  Computed  bv  the 
Axisvmmetric  and  3-D  Bubble  Codes. 
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Similar  plots  comparing  the  north  and  south  pole  displacements  in  the  two  calculations  are 
given  in  Figures  6  and  7,  respectively.  Theoretically,  the  displacements  at  bubble  maximum  in 
Figures  6  and  7  should  be  close  to  ±1 .0  when  using  nondimensional  variables.  However,  due 
to  small  inaccuracies  in  the  method’s  ability  to  predict  the  bubble’s  volume  in  these 
calculations,  the  bubble  slightly  overexpands.  This  was  expected  due  to  the  increase  in  the 
bubble’s  pressure  which  is  affected  by  the  volume  calculation. 

Figure  8  shows  a  comparison  between  the  bubble  shape  predicted  by  first-order  3-D  code 
and  the  axisymmetric  code  for  the  same  free-field  bubble  dynamics  simulations.  In  the  figure, 
the  3-D  code  shows  good  agreement  in  predicting  the  broad,  smooth  bubble  jet  typical  of 
gravity  induced  bubble  jetting.  The  3-D  code  used  a  level  four  discretization,  and  the 
axisymmetric  code  uses  nine  segments.  The  crude  discretization  applied  in  the  axisymmetric 
calculation  was  used  in  order  to  give  a  comparable  representation  of  the  bubble’s  surface 
between  the  3-D  and  the  axisymmetric  calculation. 

In  order  to  show  the  potential  for  3-D  calculations  with  the  boundary  integral  method,  a 
series  of  calculations  with  infinite  rigid  plates  at  various  angles  is  provided.  The  calculations 
include  the  effects  of  both  gravity  and  the  nearby  wall.  The  nearby  wall  is  simulated  by  use  of 
an  image  bubble.  The  figures  include  an  infinite  plate  directly  over  the  bubble  (Figure  9),  a 
plate  at  a  45°  angle  above  the  bubble  (Figure  10),  a  vertical  plate  adjacent  to  the  bubble 
(Figure  11),  a  plate  at  45°  angle  below  the  bubble  (Figure  12),  and  a  plate  directly  beneath  the 
bubble  (Figure  13).  All  of  the  calculations  performed  were  at  a  standoff  of  one  bubble  radius 
for  the  infinite  wall  and  all  calculations  used  the  same  Froude  number.  The  maximum  bubble 
radius  was  estimated  from  a  free  field  calculation  and  is  given  by  Equation  17.  The  Froude 
number  is  a  nondimensional  ratio  of  inertia  to  gravity  forces  Fr  =  V2/gL.  These  calculations 
are  for  a  Froude  number  of  8.6.  The  calculation  of  the  Froude  number  was  based  on 
Fr  =  V2/gL  =  L,/gTt2.  Additionally,  the  initial  conditions  duplicate  a  1  -lb  TNT  explosive  charge 
with  a  total  hydrostatic  head  of  33  ft.  In  these  calculations,  the  effects  of  gravity  are  fairly 
strong.  All  the  figures  are  made  using  level  four  discretization.  In  Figure  9,  the  effects  of 
gravity  and  attraction  to  the  wall  act  in  the  same  direction.  As  expected,  the  jet  is  attracted 
toward  the  wall  and  is  somewhat  sharper  and  has  a  higher  velocity  than  of  the  free  field 
calculation  (see  Figure  8),  where  the  effects  of  gravity  are  isolated.  It  is  important  to  note  that 
this  calculation  could  have  been  made  with  the  axisymmetric  code.  However,  the  calculations 
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IX 


North  Pel* 


NON-OMENSiOfMUZED  HUE 


Figure  7.  Comparison  of  Bubble  South  Pole  Surface  Displacements  as  Computed  bv  the 
Axisvmmetric  and  3-D  Bubble  Codes. 
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Axisymmetric  Calculation 


AT=1.0 


Three-Dimensional  Calculation 


lure  8.  Predicted  Bubble  Shapes  as  Computed  by  the  Axisymmetric  and  the 
3-D  Bubble  Codes. 
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Horizontal  Rigid  Plate  Above  Bubble 
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"igure  9. 


The  Calculation  of  an  Explosion  Bubble’s  Collapse  Near  an  Infinite  Rigid  Boundary 
Using  the  3-D  Bubble  Code. 


45°  Inclined  Rigid  Plate  Above  Bubble 


ATs1.95  AT-202 


Non-DImenslonalized  Time 


gure  10.  The  Calculation  of  an  Explosion  Bubble's  Collapse  Near  an  Infinite  Rigid  Boundary 
Using  the  3-D  Bubble  Code. 
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Vertical  Rigid  Plate  Along  Side  Bubble 


ATsl.O  Non-DImenslonalized  Time 


AT=2.25 


Figure  1 1 .  The  Calculation  of  an  Explosion  Bubble's  Collapse  Near  an  Infinite  Rigid  Boundary 
Using  the  3-D  Bubble  Code. 


45°  Inclined  Rigid  Plate  Below  Bubble 


Non*Dimensionaiized  Time 


Figure  12.  The  Calculation  of  an  Explosion  Bubble's  Collapse  Near  an  Infinite  Rigid  Boundary 
Using  the  3-D  Bubble  Code. 
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Horizontal  Rigid  Plate  Below  Bubble 


mm/mm/mam  mmmmm  v/mxmm  mmm. 
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Non-Dimenslonalized  Tims 


Figure  13.  The  Calculation  of  an  Explosion  Bubble's  Collapse  Near  an  Infinite  Rigid  Boundary 
Usino  the  3-D  Bubble  Codes. 

with  the  inclined  and  vertical  plate  are  not  axisymmetric  and  in  the  interest  of  consistency;  all 
of  the  figures  presented  in  this  section  are  made  using  a  full  3-D  calculation.  Figure  10  shows 
the  influence  of  a  45°  inclined  plate  above  the  bubble.  In  this  simulation,  the  influence  of  the 
plate  on  the  late  stages  of  bubble  jet  development  are  evident.  The  simulations  are  sufficient 
for  the  determination  of  jet  direction  and  magnitude  of  its  velocity.  Therefore,  this  method  can 
be  said  to  be  well  suited  for  the  study  of  bubble  jetting  in  a  strong  gravity  field  in  the  presence 
of  nearby  objects.  This  is  of  importance  from  a  military  standpoint  for  obvious  reasons.  In 
Figures  11  and  12,  the  declining  influence  of  the  infinite  plate  on  bubble  jet  directionality  can 
be  clearly  seen.  Finally,  in  Figure  13,  gravity  and  the  influence  of  the  rigid  plate  oppose  one 
another.  The  calculation  indicates  the  formation  of  a  ring  jet  around  the  lower  circumference 
of  the  bubble.  Similar  calculations  with  bubbles  between  two  infinite  plates,  where  the 
attractive  force  of  the  plates  oppose  one  another,  have  been  made  in  axyismmetric  conditions 
by  Kucera  and  Blake  (1989),  with  a  3-D  code  by  Chahine  (1982),  and  was  observed  by 
Chahine  (1989).  In  these  calculations,  the  bubble  forms  a  ring  jet  about  its  circumference 
which  eventually  separates  the  bubble,  in  an  hourglass  fashion,  into  two  bubbles.  Both 
calculations  are  similar  insofar  that  they  are  axisymmetric  in  nature  and  the  two  forces 
influencing  bubble  collapse  oppose  each  other. 
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4.  SUMMARY 


The  boundary  integral  method  was  used  to  develop  axisymmetric  and  3-D  numerical 
solutions  for  explosion  bubble  dynamics  studies.  These  formulations  were  adapted 
specifically  for  studies  involving  underwater  explosion  bubbles.  A  number  of  bubble  codes 
have  been  developed  and  are  discussed  during  the  current  formulation.  The  validity  and 
some  of  the  limitations  were  established  through  comparisons  with  previously  developed 
codes. 

The  current  effort  offers  a  number  of  alternatives  for  the  calculation  of  bubble  dynamics. 
The  3-D  code  was  unique  in  its  numerical  solution  for  the  integral  evaluation  (Wilkerson  1990). 
Further,  the  method  developed  used  a  consistent  finite  element  solution  for  interpolating 
values  of  the  potential  and  the  velocities  on  each  triangular  element. 

Future  efforts  using  this  boundary  integral  approach  will  be  centered  around  understanding 
and  controlling  inherent  instabilities  in  the  solution  process  which  are  typically  found  in 
underwater  explosion  bubble  dynamics. 
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